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Abstract. We have carried out 1.25 pc resolution MHD simulations of the ISM, on 
a Cartesian grid of < {x,y) < 1 kpc size in the galactic plane and — 10 < 2: < 10 
kpc into the halo, thus being able to fully trace the time-dependent evolution of 
the galactic fountain. The simulations show that large scale gas streams emerge, 
driven by SN explosions, which are responsible for the formation and destruction 
of shocked compressed layers. The shocked gas can have densities as high as 800 
cm~^ and lifetimes up to 15 Myr. The cold gas is distributed into filaments which 
tend to show a preferred orientation due to the anisotropy of the flow induced by 
the galactic magnetic field. Ram pressure dominates the flow in the unstable branch 
10^ <T< 10^-^ K, while for T< 100 K (stable branch) magnetic pressure takes over. 
Near supernovae thermal and ram pressures determine the dynamics of the flow. 
Up to 80% of the mass in the disk is concentrated in the thermally unstable regime 
10^ <T< 10^-^ K with ~ 30% of the disk mass enclosed in the T< 10^ K gas. The 
hot gas in contrast is controlled by the thermal pressure, since magnetic fleld lines 
are swept towards the dense compressed walls. 
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1. Introduction 

Modelling the ISM with high spatial resolution allows us to tackle a set 
of problems simultaneously, encompassing both large and small scales, 
provided that the appropriate grid size, resolution and numerical tools 
(e.g. adaptive mesh refinement) are used. To name just the most im- 
portant ones, global modelling yields information on the formation and 
lifetimes of molecular clouds, how star-forming regions are influenced 
by large-scale flows in the ISM, and which dynamical role SNe and 
superbubbles play in triggering local and global star formation. In this 
paper we investigate the effects of the magnetic field on the dynamics 
and evolution of the cold gas in the disk, the relative importance of the 
field for the ISM "phases" , and its influence on the support of shocked 
compressed layers. 
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2. Model and Simulations 

In the current work we report on kpc-scale simulations of the ISM on 
a cartesian grid of < {x,y) < 1 kpc size in the Galactic plane and 
— 10 < z < 10 kpc into the halo, using a modified version of the 3D 
SN driven ISM model of Avillez (2000) coupled to a three-dimensional 
MHD code that uses adaptive mesh refinement (AMR) in a block- 
based structure in combination with Message Passing Interface (MPI) 
suitable for massive parallel computations. 

The modified model includes a fixed gravitational field provided by 
the stars in the disk, radiative cooUng assuming optically thin gas in 
collisional ionization equilibrium, and uniform heating due to starlight. 
The radiative cooling function is a tabulated version of that shown in 
Figure 2 of Dalgarno & McCray (1972) with an ionization fraction of 0.1 
at temperatures below 10^ K and a temperature cutoff at 10 K. Back- 
ground heating due to starlight varies with z as described in Wolfire et 
al. (1995); in the midplane, at z = 0, it is chosen to initially balance 
radiative cooling at 9000 K. With the inclusion of background heating 
the gas becomes thermally stable at T< 100 K and lO^'^ <T< lO*^-^ K. 

The prime sources of mass, momentum and energy are supernovae 
whose setup is similar to that described in Avillez (2000) with some 
modifications. In the present work SNe type la with a scale height of 
325 pc are also included. In regions with a density and temperature 
thresholds of 10 cm^'^ and 100 K, respectively, an initial mass function 
is applied to determine the number of OB stars and their masses, form- 
ing an OB association. In agreement with observations 60 % of these 
stars explode within the association, while the remaining lowest mass 
stars occur in the field. The location of the field stars is determined 
kinematically by attributing to each star a random direction and a 
velocity. The associations are allowed to form in a layer with a scale 
height of 46 pc, while the field stars occur in a layer with a scale height 
of 90 pc. The time interval between the explosions of all OB stars is 
determined by their main sequence life time. 

The interstellar gas is initially setup with a density stratification 
distribution that includes the cold, cool, warm, ionized and hot gas 
"phases" in the Galaxy as described by Ferriere (1998). In these simu- 
lations the Galactic SN rate has been used and the canonical energy of 
explosion is 10^^ erg for all types of SNe. At the beginning of the simu- 
lations the uniform field components along the three axes are given by 
Bu = {Bu^o{n{z)/no)^^'^ ,0,0), where Bu,o = 3 /xG is the field strength 
and n{z) is the number density of the g function of distance 

from the Galactic midplane and uq = 1 cm~^ is the average midplane 
density. The random field component is set to zero in the beginning of 
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the simulations. This component is built up during the first millions of 
years of evolution as a result of turbulent motions, mainly induced by 
SN explosions. Self-gravity is not included in these calculations. 

The computational grid has a resolution of 10 pc, except in the layer 
between -500 and 500 pc, where three levels of AMR are used, yielding 
a finest level resolution of 1.25 pc. Periodic boundary conditions are 
applied along the four vertical boundary faces, while outflow boundary 
conditions are imposed at the top {z = 10 kpc) and bottom {z = 
—10 kpc) boundaries. All the simulations were evolved for 400 Myr 
corresponding to several sound crossing times along the fountain as 
required for the establishment of the duty disk-halo-disk cycle and the 
global dynamical equilibrium. 

3. Results 

3.1. Global Evolution 

The initial evolution of the magnetized disk is similar to that seen in 
the HD runs (Avillez 2000, Avillez & Breitschwerdt 2003), that is, the 
initial stratified distribution does not hold for long as a result of the lack 
of equilibrium between gravity and (thermal, kinetic and turbulent) 
pressure during the "switch-on phase" of SN activity. As a consequence 
the gas in the upper and bottom parts of the grid collapses into the 
midplanc, leaving low density material in its place. However, in the 
MHD run it takes a longer time for the collapse to be completed as 
a result of the magnetic pressure and tension forces. As soon as the 
system has collapsed and enough supernovae have gone off in the disk 
building up the required pressure support, transport into the halo is 
not prevented, although the escape of the gas takes a few tens of Myr 
to occur. The crucial point is that a huge thermal overpressure due 
to combined supernova explosions can sweep the magnetic field into 
dense filaments and punch holes into the extended warm and ionized 
Hi layers. Once such pressure release valves have been set up, there 
is no way from keeping the hot over-pressured plasma to follow the 
pressure gradient into the halo. As a consequence the duty disk-halo- 
disk cycle of the hot gas is fully established, which combined to the 
input of energy into the ISM by SNe, diffuse heating, energy lost by 
radiative cooling and magnetic pressure leads the system to evolve into 
a dynamical equilibrium state within some 200 Myr (Kahn 1981). 

Fig. 1 shows slices of the three-dimensional data cube of the density 
and magnetic field distributions in the Galactic midplane. The highest 
density (n > 100 cm~^) gas tends to be confined into shocked com- 
pressed layers that form in regions where several large scale streams of 
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Left panel: mavinez_figla, Right panel: mavillez_figlb 



Figure 1. Density and magnetic field distribution in the Galactic midplane after 374 
Myr of disk evolution. The resolution of the finest AMR level is 1.25 pc. 

convergent flow (associated with laminar flows driven by SNe) occur. 
The compressed regions are fllamentary in structure, tend to be aligned 
with the local field and are associated with the highest field strengths^. 
The formation time of these high density structures depends on how 
much mass is carried by the convergent flows, how strong compression 
is and on the rate of cooling of the compressed region. The compressed 
layers have on average lifetimes of 10-15 Myr. The streams that pro- 
mote the formation of the shocked dense layers are also responsible 
for their destruction. The filamentary structures are thicker in these 
simulations than in similar HD simulations, because the compression 
of a shock wave at given Mach number must go into increase of thermal 
and magnetic pressure in the MHD case and only into thermal in HD. 

3.2. Field Dependence with Density 

The simulations started with a fixed value for the magnetic field varying 
with z and having a midplane value of 3 x 10~® G. During the evolution 
of the system thermal and dynamical processes broaden the distribution 
of the field strength in such a way that after the global dynamical 
equilibrium has been set up the field strength in the disk spans two 
orders of magnitude from 10~^ to 10~^ G as can be seen in Fig. 2. The 
figure shows a scatter plot of the field strength as function of density 
in the simulated disk after 400 Myr of evolution. 

^ A similar result is seen in similar HD simulations. However, there is no preferable 
orientation for the filaments. 
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Left panel: mavillez_fig2a, Right panel: mavillez_fig2b 



Figure 2. Scatter plot of the magnetic field strength as a function of density for the 
T< 10^ (red) and 10* <T< lO^^'^ K (black) regimes (left panel) and 10^ <T< 10* 
(black) and T> 10®-^ K (red) regimes (right panel) at 400 Myr of evolution. 

Both panels in the figure show a large scatter in the field for the 
same density. This lack of correlation suggests that the field may not 
follow the "classical" relation of B ^ p^/^. This essentially means that 
the equation of state is considerably softer than in the pure magnetic 
case and hence the effective adiabatic index 7 must be well below 
unity. This may be attributed to the variety of processes that affect the 
energy balance of the magnetized plasma, such as turbulence, heating 
and cooling. Specifically the equation of state becomes stiffer, if the 
dependence of the heating rate on density is steeper than that of the 
cooling rate. 

The large scatter in the field seen in Fig. 2 also suggests that the 
field is being driven by the inertial motions, rather than it being the 
agent determining the motions. In the latter case the field would not 
be strongly distorted, and it would direct the motions predominantly 
along the field lines. The high field variability is also seen in the right 
panel of Fig. 1, which shows a highly turbulent field, that seems to be 
uncorrelated with the density. 

3.3. Pressure Variation with Temperature 

Left panel of Fig. 3 shows the scatter points of the thermal (black) 

and magnetic (red) pressures as a function of temperature, while the 
right panel shows the average magnetic, thermal and ram pressures as 
function of temperature. From the figure it is seen that all the T< 10^ K 
gas has Pb > Pram ^ Pth : demonstrating that magnetically dominated 
regions do exist, while the hot gas has the highest thermal pressure and 
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Left panel: mavillez_fig3a, Right panel: mavillez_fig3b 



Figure 3. Loft panel shows a scatter plot of the magnetic (rod) and thermal (black) 
pressures as functions of temperature at 400 Myr. Right panel shows the average 
magnetic (red), thermal (black) and ram (green) pressures as functions of temper- 
ature at 400 Myr. This plot is obtained by calculating the average pressure in a 
certain temperature bin. The red and black curves represent the average values of 
the pressures seen in the left panel. 

the lowest magnetic pressure. For T< 10^ the relation B cc does not 
hold, as a varies between -0.06 and 0.085 (Avillez & Breitschwerdt 
2004) indicating that the magnetic and thermal pressures are inde- 
pendent. A result that is consistent with the almost zero variation of 
magnetic pressure with temperature seen in the right panel of Fig. 3. 

The figure also shows that for 10^ <T< 10^ K ram pressure domi- 
nates over magnetic and thermal pressures, that is, ram pressure deter- 
mines the dynamics of the flow, and therefore, the magnetic pressure 
does not act as a restoring force as it was already suggested by the lack 
of correlation between the field strength and the density (for a detailed 
discussion see Passot &; Vazqucz-Semadeni 2003). The alternating dom- 
inance between magnetic and ram pressures seen for 10^ <T< 3 x 10^ 
K results from the fact that about half of the explosion energy in 
supernovae becomes kinetic, and therefore, ram pressure becomes also 
important near the energy sources. 

The dynamical picture that emerges from these simulations is that 
thermal pressure gradients dominate mostly in the neighborhood of 
supernovae, which drive motions whose ram pressures are dominant 
over the mean thermal pressure (away from the energy sources) and the 
magnetic pressure. The magnetic field is dynamically important at low 
temperatures. The field does not act only by its pressure gradient, but it 
acts dynamically also by magnetic tension forces, l/47r(BV)-B. Whereas 
magnetic pressure is isotropic, magnetic tension is always along the 
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Left panel: mavillez_fig4a, Right panel: mavillez_fig4b 



Figure 4- Averaged mass weighted (left panel) and cumulative (right panel) his- 
tograms of the density in the Galactic disk calculated between 300 and 400 Myr 
using 101 snapshots with 1 Myr of interval. The black and blue solid lines represent 
gas in the thermally stable regimes at T< 10^ K and 10^ '' <T< lO"'^ K, the orange 
solid and dashed lines represent the thermally unstable gas with 10^ <T< 10^ K 
and 10^ <T< lO^ '' K, respectively. The gas with lO"*'^ <T< lO^ '^ K is represented 
by the green line, while the hot (T> 10^'^ K) gas is represented by the red line. 

lines of force. Equilibrium only exists in force-free configurations for 
which V X B = 0. However, if there is a thermal pressure then it is 

practically impossible to balance the magnetic stresses and most of the 
configurations will be unstable. In other words, in a dynamical system 
like the ISM, which includes systematic gas motions, j x S-forces are 
inevitable, and hence the magnetic field can be dominant. Thus, regions 
with a high field strength such as the cold clouds should be magnetically 
controlled. 

3.4. Distribution of the ISM Mass 

Up to 80% of the ISM mass in the simulated disk is concentrated in 
the unstable regime at 10^ <T< 10^-^ K (orange curves in Fig. 4) 
where '-^ 30% of the ISM mass is concentrated in 10^ <T< 10^ K gas 
(dashed orange curve), while the 10^ <T< lO^-S K gas (sohd orange 
line) encloses some 50% of the mass. The stable regimes T< 10^ K 
and lO^ '^ <T< lO^ ^ K enclose ~ 1% and ~ 10% of the disk mass, 
respectively. Note that most of the cold gas is concentrated in the 
unstable branch at temperatures T< 10^ K and has densities up to 
100 cm"^ while gas with highest density has T> 10^ K. The hot gas 
amounts to less than 1%, while the unstable regime with 10^"^ <T< 
10^'^ encloses ^ 10% of the total mass. 
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4. Summary 

The highest resolution MHD simulation of the ISM carried out to date, 
discussed in this paper, shows that large scale streams driven by SNe 
are responsible for the formation and destruction of high density clouds 
in shocked compressed layers. The clouds formed within the shocked 
gas can have densities as high as 800 cm~^ and lifetimes up to 15 Myr. 
Ram pressure dominates most of the coolest flows (in the unstable 
branch 10^ <T< 10^-^ K) except at temperatures below 100 K (in the 
stable branch), where magnetic pressure takes over. 

Most of the mass in the disk is concentrated in the thermally un- 
stable regime for T< 10^-^ K with some 30% of the mass enclosed in 
gas with T< 10^ K. This is consistent with observations and numerical 
simulations discussed by Heiles (2001), Gazol et al. (2001) and Kritsuk 

Norman (2002). However, these results are certainly in disagreement 
with classical ISM theories (McKee &: Ostriker 1977), as from our sim- 
ulations up to 90% of the ISM mass is found to be in the thermally 
unstable regimes with 10^ <T< lO^-^ K (~ 80%) and lO^-^ <T< lO^-^ 
K (~ 10%), where it should not exist, if it were not for the importance of 
dynamical processes. This calls into question the standard paradigm of 
an ISM distributed over three phases and being in pressure equilibrium. 
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